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Abstract 

Since computational efficiency and wave resolution scale with accuracy, the ideal would 
be infinitely high accuracy for problems with widely varying wavelength scales. Currently, 
many of the computational acroacou sties methods are limited to A Ui order accurate Runge- 
Kutta methods in time which limits their resolution and efficiency. However, a new proce- 
dure for implementing the Modified Expansion Solution Approximation (MESA) schemes, 
based upon Hcrmitian divided differences, is presented which extends the effective accuracy 
of the MESA schemes to 57 /7 ' order in space and time when using 128 bit floating point 
precision. This new approach has the advantages of reducing round-off error, being easy 
to program, and is more computationally efficient when compared to previous approaches. 
Its accuracy is limited only by the floating point hardware. The advantages of this new 
approach are demonstrated by solving the linearized Euler equations in an open bi-periodic 
domain. A 50(V 7 ' order MESA scheme can now be created in seconds, making these schemes 
ideally suited for the next generation of high performance 256-bit (double quadruple) or 
higher precision computers. This ease of creation makes it possible to adapt t he algorithm 
to the mesh in time instead of its converse: this is ideal for resolving varying wavelength 
scales which occur in noise generation simulations. And finally, the sources of round-off 
error which effect the very high order methods are examined mid remedies provided that 
effectively increase the accuracy of the MESA schemes while using current computer tech- 
nology. 


1 Introduction 

Predicting the sources of jet noise requires comput at ional methods that are orders of magnit ude 
more efficient and that provide* very high resolution. This is accomplished numerically with very 
high accuracy, adaptable, explicit methods whose benefits are as follows: 

. High accuracy methods arc* more efficient and provide finer resolution of the physics [5]; 

• Adaptable methods can adjust their accuracy to resolve steep gradients while avoiding the 
complexities of mesh adaptation; 

• And. explicit methods permit highly parallel/scalable computations by minimizing inter- 
processor communication [4]. 

The proposed approach enables the MESA schemes [5] to accomplish those objectives by mak- 
ing them simple* to program, adapt, and compile while simultaneously reducing floating point 
operations and round-off error. 

MESA schemes can be viewed as a multidimensional, higher order extension of Lax-Wendroff 
scheme's that incorporate* more of the physics (via cross-derivative* information) necessary to more* 
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accurately propagate waves along their characteristic surfaces. The MESA schemes require' es- 
sentially two procedures; a spatial interpolation followed by a time advance. Previously, the 
second step was implemented using a recursive definition [4] which enabled arbitrarily high ac- 
curacy in time. However, the spatial interpolation step required computer algebra [3] for the 
symbolic creation of one-dimensional intcrpolants. effectively limiting the accuracy of MESA 
schemes to 29 th order because of the computational complexity in producing its symbolic form. 
This final limitation is removed in this work by replacing the one-dimensional symbolic in- 
terpolant with a simple and efficient form of the Hermitian divided difference interpolant. A 
key finding of this paper is that all the spatial derivatives of a Hermitian divided difference 
interpolant. at the midpoint of a two-point, multidimensional stencil, have a simple algebraic 
expression which eliminates the need for computer algebra tools. Since small stencils have many 
advantages [7] such as better resolution and ease of boundary implementations, using this new 
form of Hermitian divided differences with two-point Hermitian MESA schemes would appear 
ideal. 

Divided differences have been used to interpolate data for many years dating back to Isaac 
Newton. [11] and [2], With the advent of digital computers, divided differences have been re- 
placed by splines since polynomial interpolations tend to oscillate at higher orders. In addition. 
Hermitian methods have not been extensively used due to the difficulty of obtaining the deriva- 
tives of the function being approximated [14]. However, in this work, polynomial oscillations 
are eliminated since only a single interval is used and the derivatives of the function at the end- 
points of this interval are completely prescribed. Also, round-off error is reduced for high-order 
interpolations since guard digits are introduced into the tableau [13]. And Hermitian divided 
difference interpolations coincidentally use the same data found in the stencil of a 2* t s 1 order 

MESA scheme (/. } * ). For these reasons. Hermitian (Birkhoff [10]) divided 

difference's are used to great advantage' here. 

This paper first describes the new approach to spatial interpolation in one-dimension and 
then extends it to 2 x 2 Hermitian stencils. Next, the linearized Euler equations are solved in 
a bi- periodic open domain by applying this new interpolation method to the MESA schemes. 
The error and efficiency of the now approach is then compared with the previously best known 
approach. And finally, various techniques are shown for improving the effective accuracy of very 
high order methods (> 30) when computer precision is limited. 


2 Two-Point One-Dimensional Hermitian Divided-Difference 
Interpolation 


We will now provide an overview of Newton's interpolation method based upon divided differ- 
ences. 

Let x, : / = 0. 1 n be any (n+l) distinct points of [a.b] and let f be a differentiable 

function. C d [a*b]. The coefficient of x n in the polynomial p e P n that satisfies the conditions 


dp(Xj) 

Ox 


Of(Xj) 

dx 


= 0, 1 //. 


( 1 ) 


is defined to be a divided difference of order n [13]. These divided differences, once deter- 
mined. completely define the spatial interpolant satisfying equation 1. 

A convenient mnemonic commonly used to determine the divided differences is shown in 
figure 3 and is referred to as a divided difference tableau. In this figure is the tableau for 
a fifth order one-dimensional interpolant. WTien Hermitian data (the data and their spatial 
derivative's) are used for interpolation [13]. the tableau is formed by inserting all the known data 
into each column of the tableau as shown in figure 4 for the case of a two-point Hermitian stencil 
with throe data elements per grid point (f(x< y), f x {x. y). f xx (x. •*/)). In general, for a two-point 


NAS A/TM— 2000-209944 





G 1 O 

X = Xu X = 0 X = Xi 


Figure* 1: Two Point One-Dimensional Stencil 


X = Xu x = 0 x = X\ 
V = 2/i V = 2/i 



X = Xq X = 0 X = Xi 

2/ = 2/o iJ ~ Vu 


Figure 2: Two Point Two-Dimensional Stencil 


Hermit. ian stencil in which tlio data: f {x. y). f x (x< y). f xx (x.y) i‘ s known at both 

grid points (x=x 0 and x— Xi) in figure 1. the following procedure will correctly place this data 
into the tableau. 

Do[Do[Q[iJ} = ^±;Q[i + « + Lj} = .{jA).i}}. {Ul,}}; (2) 

With this initial data inserted into the tableau and the distance between both points defined 
by A/i, the rest of the tableau is constructed using [2]: 



Figure 3: Divided-Difference Tableau for c2o2 MESA scheme* 
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figure 4: Loaded Divided-Difference Tableau for c2o2 MESA scheme. X-Direction 


Do[Do[Q[i.j] = 


(QjLj - 1] - Q[i - Lj - 1]) 


A/; 


. {j. i - ». <}]. {/. * 1. 2 * .s + 1 }]; 


(3) 


And with this tableau, the interpolating Hcrinitian polynomial on a Hermitian stencil with 
two points and .s • 1 data elements (primitive variable and its spatial derivatives) per grid point 
may be evaluated with Newton's interpolators’ divided difference formula [2]: 


t>«(s+l)-l ,-i 

f(x) = IJf 1 - Xj) 

i-o j - o 


(• 1 ) 


which may bo rewritten as: 


fix) = y^Q(i.i)(x - Xo)' 4- Q(i. i)(x - x 0 ) s+1 (i - Xi) ,-< ® +1) (5) 

/— 0 i—s f 1 

Hermitian MESA schemes require evaluating this polynomial at the center of each stencil 
to interpolate the solution data and their spatial derivatives. In the general case in which a 
MESA scheme of arbitrary accuracy is used, the spatial derivatives of equation 4 become very 
complicated and require computer algebra tools for their solution. This limits the accuracy and 
adaptability of these methods. 

For example, a 49 //l order MESA scheme in one-dimension requires 25 data elements per grid 
point. But to time advance all those elements requires determining the derivatives of equation 4 
up to the i)Q fh order. The product terms in that equation will double after each differentiation 
by the product rule of differentiation; this produces equations with approximately 2 50 % 10 1(i 
terms. Besides quickly exhausting memory resources of most computers, this also takes too 
much time to calculate' and makes it difficult to modify the accuracy of the numerical scheme 
in real time to accommodate steep gradients. 

Fortunately, by using a two-point stencil and a reformulation of equation 5. it is possible to 
efficiently calculate all the necessary spatial derivatives without the use of computer algebra. 


3 Fundamental Result: Direct Interpolation of Spatial Deriva- 
tives at Center of Two Point Hermitian Stencil 

A 2*s- 1 order Hermitian MESA scheme (labeled c2os) will contain (,s~ !).(* - I) 2 , or (,s + ] ) 3 data 

elements per grid point for each primitive variable in one. two. or throe-dimensions respectively. 

The MESA scheme requires all these data (‘lenient s to be advanced in time. Accomplishing this 
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requires interpolating 2(s4- l),4(.s + l) 2 , or 8(.s+ l) 3 spat ial derivatives at the center of the stencil 
for each primitive variable in one. two, or three spatial dimensions respectively. Each of those 
interpolations normally requires evaluating the derivative of equation 4 at x = 0. Because of the 
product term in this equation, higher order derivatives of this equation become' complicated, as 
mentioned. However, the following main result of this paper provides an alternative* formulation 
for the interpolation of the spatial derivatives at the center of a two-point Hermit ian stencil. 


d dx f(x = 0 ) 
dx dx 


£ 


Q(i-i) 


i—dx 

where Z is defined as: 

dx r 


(/ — dx)\ 


, r ,\(*-dx) 


( — J'o) 


2(*+l )-J 


y [Q(i' i)Z{L s. dx, x 0 . Xi )] (6) 


*=#+ 1 


Z(i. .s. dx. Xq, Xi ) = ^ 

r— 0 

with P[ and P> defined as: 


dxl 


{dx — r)!r! 


(_ Io ) ( »+i-' )(_ Xl )('-(^ n-d*+,-) Pi ,.)p 2 (. s . r) 


dx—l—r 


Pi(i. #, dx. r) = JJ [* - (» ~ 1) - 1 


and 


7" — 1 


1-A-] 


7) 


( 8 ) 


(9) 


.respectively. 

This form can be rewritten as : 


k—0 


2 s+1 




( 10 ) 


i—dx 


where the function coef(i.dx) is predefined as: 

Do[Do[corf[i.dx) = ( , - jb ^ T (-x 0 ) ( '- dj - ) . (i. ( /x..s}]. {(/x. 0. (2* + 1)}] 

Do[Do[ 

coc f[Ldx\ = 

Eto[(rf^WT(-^o) l * +1 -' ) (-ii) ( '-*- 1 -^ + '' ) * (ID 

Product[(i — ,s — 1 — e). {7 , 0. dx — 1 — r }] * 

Product [(.s -f 1 — A*), {A*. 0. r — 1 }]] 

1, 2s -rl}]. {dx. 0. (2s + 1)}]; 

The function coef(i.dx) in equation 10 is independent of spare and time and needs computed 
only once. Thereafter only the divided difference terms Q(i,i) must be evaluated for each stencil 
at each time step using equation 3. 


4 Cost Comparison 

As mentioned, this new approach eliminates the need for using computer algebra to produce 
high order MESA schemes; this reduces the initial cost of programming the MESA schemes and 
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increase's its possible accuracy and adaptability. In addition, another unanticipated benefit of 
this new approach is fewer floating point operations are required as compared to the best previous 
approaches. For example, the total cost to evaluate all of the spatial derivatives required in a 
2 * .s + 1 order MESA scheme (c2os) is the sum of the cost to compute the divided difference 

terms, Qs, using equation 3 and the cost of evaluating all of the spatial derivatives ( ° ) 

using equation 10. Therefore the total cost of the new approach is: 

((1 + «s ) 2 ) + (3 5.s- -t 2s 2 ) = 4 -r 7s — 3s 2 (12) 


multiplications. 

In the previously best known procedure [4], in which computer algebra is used to generate' 
the one-dimensional interpolant. each spatial derivative required at most 2s -t- 2 multiplications, 
and there were 2.s + 2 terms to bo evaluated in one dimension resulting in a total cost of: 

(2* -r 2) 2 = 4 + 8# + 4* 2 (13) 


multiplications. 

Therefore the new Hermit ian (Birkhoff) divided difference form requires 

a 4 .s 2 (14) 

fewer multiplications. However, the assumption that calculating each spatial derivative under 
the old approach required 2s - 2 operations is an upper limit and in practice certain algebraic 
cancellations may reduce that number. At higher accuracy this upper limit is more likely since 
the equations become' large and cancellations are more difficult. 


5 Two-Point Two-Dimensional Hermitian Divided-Difference 
Interpolation 

It is necessary to perform a multidimensional interpolation for the MESA schemes in two and 
three dimensions. While this could possibly be accomplished using other multidimensional 
divided difference techniques [12]. we will use the tensor product approach as in Dyson [4] since 
this permits a detailed comparison of the new and old approaches. 


5.1 Tensor Product Approach Overview 

The tensor product approach interpolates all the spatial derivatives required for the MESA 
schemes by performing a series of one-dimensional interpolations [4]. Each one-dimensional 
interpolation requires a new divided difference tableau to be generated using equation 3. A set. 
of one-dimensional interpolations is first performed in the x-direction to interpolate the data 

- — : i = 0. 1 s. and j = 0. 1.2 2 * « -r 1 at y = yo and y — y\. And then. 

by interpolating in the y-direction using only these interpolated values, the following terms are 
found: 


d i +Jf(T = [ly = ()) 
dx'yj 


VLj:iJ = 0.1,2 2*.s+ 1 


(15) 


using the coordinate system in figure 2. These terms are required for time advancing all the 
data on each grid point, namely: 


f(x.y) 

dx'yj 


VLj 


Lj = 0 . 1,2 * 


(16) 
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at grid points (x 0 . yo)> (#o* Vi h (*i ■ ?/u)< ftnd(xi. iji) as shown in figure 2. 

This tensor product procedure is extensible to any size stencil, but the Hermit ian divided 
difference one-dimensional interpolat ion as developed in this paper is limited to two-point inter- 
polations. 


5.2 Fifth Order Two-Dimensional Interpolation Example 

The combination of one-dimensional two-point Hcrmitian divided difference interpolation with 
multidimensional tensor product extensions is best understood with a simple example. We will 
completely describe the process required to interpolate the data on a two point two-dimensional 
stencil shown in figure 2 for the 3 th order MESA scheme. c2o2. 

The c.2o2 MESA scheme in two dimensions, will contain the data: /. /. r - /.r*- fy fyx - fyxx * 
fyy. f yyx , fyy XX at all four grid point locations as shown in figure 2 for each primitive variable 
(pressure, u- velocity, v- velocity). Using t hose 3G pieces of information for each primitive variable', 
the two-dimensional spatial interpolant of the form: 

5 o 

f(x.y) = ^2^2rf{Lj)x’y J (17) 

7-0 j = 0 


is created by finding the values of the cf(Lj) coefficients where 


rf(Lj) 


fdJ) 


(18) 


These coefficients are found by first performing one-dimensional interpolations of the form: 


fix) 5>/(<V (19) 

7—0 


where 


cfd) = 


f {i) (x = iky) _ df(x = 0. y ) ( 1 

i\ Ox 1 


( 20 ) 


By letting f(x) in equation 19 be replaced by /(r). f y (x ). and f yy (x ). respectively, then 

= 0 .y) 1 0 {nj) f(x = 0) ( 1 


<‘fd) = 


Oyj 


dx l 0y J 


( 21 ) 


for j=0. 1. and 2. respectively; And. i=0,l,2.3.4. and 5. 

Thus, it is possible to interpolate all the values, cf(i) for j=0.1. and 2 at the location 
[x = 0. y = y 0 ) shown in figure 2 using the three one-dimensional interpolants of the form 
of equation 19. This is repeated for the cf(i) values at the location (r = (J. y = i/i) for each j=(). 
1. and 2 for a total of six one-dimensional interpolations. 

At this point, we have the two sets of interpolated values: 


d*+ J f(x = ny) n\ 

0x } dy J V / ! / 


( 22 ) 


for y — ijo and y i; and for i=0.1.2,3.4.5; and j=0.1, and 2. 

This data may be used to evaluate the derivatives in equation 15 at the center of the stencil. 
This requires a second set of one-dimensional divided difference interpolations in the v-direction 
using only this data. 
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Interpolating in the y-dircction is accomplished with the polynomial: 


where 


' 70 ') = 


f{x = 0. y) = Y^('fU)y J 
j = 0 

d J f(x = 0 .y = 0) ( 1 


( 23 ) 


(24) 


dy J \j- 

This time however, we will interpolate the functions: f -fx-fxx-fxxx-fxxxx- and fxxxxx by 
performing the following substitution into equation 24: 


fix = O.i/) = 


d’fix = 0 .y) ( 1 


dx’ 


i) 


for i=0. 1 .2.3.4. and 5 so that 


rf(j) = 


d'+ J fix = O. y = 0) f 1 


dx l dyj 


i\j\ 


(25) 


(26) 


Therefore this requires six more one-dimensional divided difference interpolations for a total 
of twelve when including the six horizontal interpolations. The substitution of equation 25 into 
equation 24 reuses the data from the x-direction interpolation, thereby introducing efficiencies 
not found with other multidimensional interpolation procedures. Note that this efficiency is 
only possible if the data in equation 16 are available at each grid point so that a symmetry of 
the spatial derivatives exists in all dimensions because both the tensor product and Hermitian 
divided differences require this complete set of derivative information. Fortunately, the two-point 
Hermitian MESA schemes (c2os) provide exactly this information. 


5.2.1 Horizontal Interpolation Procedures 

Applying the above concepts can be reduced to a few simple steps. Interpolating the intermediate 
values at (0. ijo) and (0. tq) in figure 2 is accomplished by: 

• First, loading the known data from the stencil into the tableau shown in figure 4 to 
interpolate the data at y = y o. The tableau is loaded with the function /. f x . and 
data contained at the two grid point locations x 0 and jq as indicated in figure 1. 

• Second, build the divided difference tableau as described in equation 3. 

• Third, evaluate the spatial derivatives '* f or i— 0.1 .2.3.4. and 5 using equa- 

tion 10. 


• Fourth, repeat these three steps by substituting f(x) with / = /. / y . and f uu . 

• Fifth, repeat these four steps at y = y \ . 

After these procedures, wo have calculated the data: 

O iiJ fix = 0 ,y) w . . , 0 . , w • • ,> , „ 

— — Vt : i = 1 . 2 . 3 . 4 . otind V j : j = 0 . 1. 2 

dx'dyJ 

at grid coordinates (0. ij u ) and (0. t/i ) as labeled in figure 2. 


(27) 
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Figure 5: Loaded Divided-Difference Tableau for c2o2 MESA scheme. Y- Direct ion 

5.2.2 Vertical Interpolation Procedures 

Next, the data shown in equation 27 is used to perform the Henriitian divided difference' inter- 
polation in the v- direct ion by: 

• First loading that data into the tableau as shown in figure 5 to interpolate the intermediate 
data of the last step along the line x = 0. The tableau is loaded with the function /. 
and f yy data contained at the two grid point locations (0, ?/o) and (O.jf/i) as indicated in 
figure 2. 

• Second, build the divided difference tableau as described in equation 3. 

• Third, evaluate the spatial derivatives for j=0.1.2,3.4, and 5 using equation 10. 

• Fourth, repeat these three steps by substituting f{x) with / — f ' I x ■ f XX - fxXX' fx, XXX' 
and fxxx xx which are the intermediate data previously evaluated using the horizontal 
interpolation procedures. 

After these steps, all the spatial derivatives of f(x,y) necessary for advancing the data on the 
grid is available for the MESA scheme. 


6 Results 


For the purposes of comparing this new approach with the old approach, it is necessary to divide 
the results from the new approach by a factorial term. This is because the coefficients for the 
previous method [4] were equivalent to using the equation: 

2*.s+l 

/(x) = ^2 cf{i)x' (28) 

/'=() 


in which the factorial term is included in those cf(i) coefficients instead of the form: 

2 * S -j- 1 r / • \ 

<■/( > ) ^ 

/-() 


m = v 


(29) 


in which the factorial term is not included and therefore cf(i) = f) .{£? ■ . 

W ith t his correction factor a direct comparison is possible 1 between both approaches for very 
high accuracy by using automatic code generation [6]. In particular, the relative efficiency and 
accuracy of applying both interpolation methods to the same wave* propagation problem are 
measured. 
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6.1 Problem Definition 


Wave propagation is described by the linearized Euler equations and its correct simulation in 
time is important in many aeroacoustic applications. We will solve the bi-periodic open domain 
problem in which the physical domain is a unit square ([-1. 1] x [-1. 1] x [0. 3~]). The solution of 
the linearized Euler equations in this case is assumed to be y-periodic (top and bottom of box 
repeat) and x-periodie (left and right sides of box repeat). Using separation of variables with 
periodic boundary conditions, on the following linearized Euler equation system with a constant 
mean convection velocity vector ( M x .M y ): 


dp 

dt 


du 

dt 

dv 

dt 

M x 


. , du , , du dp 

M x — + M y — -r = 0. 

dx dy dx 


M, 


dv 

dx 


dp 

dx 


, , dv dp 

Iu dy ~dy = 

+ v - v = '»• 

dy dx dy 


(30) 


with the boundary conditions : 


p(l.y.t) =]>{-\.y.t) 
u(l.y.t) = u(-l.y.t) 
v(l.y.t) = v(-l.y.t) 
p(x. 1.0 = p(x. -1.0 
u{x. 1.0 = n(x . —1.0 
v(x. 1. 0 = «’(x. —1.0 


provides the following analytical solution: 

p(x.y.t) = cos (~tV2) sin(^-(- (M x t) x)) sin(~(- {M y t) »- y)) 
u{x f) = cos(x (— (M x t) • x)) sin(/tt\/2) sin(7r (— (fl/^Q — y)] 

V2 

r(x y 0 = cos (~ (Myt) ~ </)) sin(7rty/2) sinpr (- (Al x t) + x)) 

y/2 


(31) 

(32) 

(33) 


6.2 Numerical Results 

The results of these numerical experiments with 8 grid points per wavelength are shown for 
double- precision (64- bit reals) in tables 1 and 2 and for quadruple precision (128-bit reals) in 
table 8. The time to complete* the simulation using each approach is a measure of relative effi- 
ciency and is shown for each MESA scheme tested. In addition, the results for both approaches 
arc* placed above and below each other for easy comparison. The seven column headings are 
defined as: 

N Number of time steps. 

T Total time elapsed in simulation. 

niaxperr The maximum absolute error in the pressure. 

Uperr A measure of the average error (Ll norm) . 
phmax Maximum pressure occurring in domain. 
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phniin Minimum pressure occurring in domain. 

e-ratio The change in the total energy content of the system as a ratio 7 ^ 77^77 (should be one 
for no change). 

Despite the cost analysis in section 4 showing a slight advantage to the new approach for 
all MESA schemes, numerical experiments show the old approach is actually faster for methods 
less than about 15 f/l order. This is likely due to the length of the one-dimensional equations 
being less than the worst case used in the analysis. However, for higher order methods the new 
approach is modestly faster while using double-precision calculations and it is significantly faster 
when quadruple-precision is used (see table 3). 

Notice in tabic' 1 that, the new approach begins demonstrating less round-off error at, ap- 
proximately ll " 1 order accuracy. And at 15 th order accuracy the new approach still maintains 
essentially no growth in the energy compared to the old approach in table 2. And at \7 th order 
both methods begin introducing significant round-off errors into the overall energy. 

In quadruple precision, both approaches appear similar at \o th order, however the new 
approach is about 50 percent faster. And at 21*' order accuracy the round-off errors are higher 
in the old approach while the new approach is more efficient. By the time 27 th order accuracy 
is reached, both approaches are showing significant round-off error in the total system energy. 

It was not possible to create a 33 rd order or higher MESA scheme using the old approach, 
but the new approach appears to maintain accuracy up to 57' /j order in quadruple' precision 
before round-off error becomes too excessive. The sources of this round-off error are explained 
in the appendix A. 
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N T 

maxperr 

llpcrr 

phmax 

phmin 

e-rat io 

OLD APPROACH ABOVE. NEW APPROACH BELOW 

c2o2: Old Approach Timc=254.79 , New Approach Time=290.39 

20U 

200 

1.0U000E+01 

l.UOOOUE+Ol 

9.07553E-04 

9.07553E-04 

1.36338E-03 

1.36338E-03 

8.95490E-01 

8.95490E-01 

-8.95490E-01 
-8.954 90E-01 

9.97694E-01 

9.97694E-01 

ffwllB 


2.06528E-03 

2.06528E-03 

3.24237E-03 

3.24237E-03 

2.40990E-01 

2.40990E-01 

-2.40990E-01 
-2.40990 E-01 

9.79382E-01 

9.79382E-01 

LO £0 
11 

1.00000E+03 

1.00000E+03 

7.401 71 E-02 
7.401 71 E-02 

1.1 1612E-01 
1.1 1612E-01 

7.04408E-01 

7.04408E-01 

-7.04408E-01 

-7.04408E-01 

8. 1406 IE-01 
8. 1406 IE- 01 

c2o3: Old Approach Timc=572.39, New Approach Timc=641.52 

200 

200 

1.00000E-01 

1.00000E+01 

2.04781E-06 
2. 0478 IE- 06 

2.99855E-06 

2.99855E-06 

8.96396E-01 

8.96396E-01 

-8.96396E-01 

-8.96396E-01 

9.99996E-01 

9.99996E-01 

2000 

2000 

1.00000E+02 

1.00000E+02 

5.51277E-06 

5.51277E-06 

1.08654E-05 

1.08654E-05 

2.43050E-01 

2.43050E-01 

-2.43050E-01 

-2.43050E-01 

9.99932E-01 

9.99952E-01 

20000 

20000 

1.0U000E-t03 

1.00000E-03 

1.77378E-04 

1.77378E-04 

2.67937E-04 

2.67937E-04 

7.78247E-01 

7.78247E-01 

-7.78247E-01 

-7.78247E-01 

9.99543E-01 

9.99543E-01 

c2o4: Old Approach Time— 1108.63, New Approach Tinier 131 1.10 | 

HI 

1.00000E+01 

l.GOOOOE-Ol 

2.86266E-09 

2.86266E-09 

4.24757E-09 
4.2475 7E-09 

8.96398E-01 

8.96398E-01 

-8.96398E-01 

-8.96398E-01 

I.OOOOOEtOO 

1.00000E+00 



8.95464E-09 

8.95463E-09 

1.86641E-08 

1.86641E-08 

2.43055E-01 

2.43055E-01 

-2.43055E-01 

-2.43055E-01 

1.00000E-00 

1.00000E+00 



2.45441E-07 
2.4544 IE-07 

3.841 70E-07 
3.84169E-07 

7.78425E-01 

7.78425E-01 

-7.78425E-01 

-7.78425E-01 

9.99999E-01 

9.99999E-01 

c2oo: Old Approach Tiine=2034.22. New Approach Time=2533.79 

m 


2.70040E-12 

2.69773E-12 

4.09499E-12 

4.09191E-12 

8.96398E-01 

8.96398E-01 

-8.96398E-01 

-8.96398E-01 

9.99999E-01 
1.00000E— 00 


1.00000E-02 
1.00000E -02 

2.00462E-1 1 
2.00432E- 1 1 

3.454 19E- 11 
3.45456E- 1 1 

2.43055E-01 

2.43055E-01 

-2.43055E-01 

-2.43055E-01 

1.00000E+00 
1 -00000E+00 

20000 

20000 

l.OUOOOE - 03 
1.00000E— 03 

1.53291E-09 

1.53274E-09 

2.87923E-09 

2.87917E-09 

7.78425E-01 

7.78425E-01 

-7.78425E-01 

-7.78425E-01 

1 .00000E • 00 
I.OOOOOETOO 


Tabic 1: Cost and Round-off Error Comparison of New and Old 2D Interpolation Methods With 
8 Grid Points Pei Wavelength: Double Precision Computations 
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N 

T 

maxperr 

llperr 

phiriax 

phmin 

e-ratio 

OLD APPROACH ABOVE. NEW APPROACH BELOW | 

c2o6: Old Approach Time=3372.85, New Approach Time=3978.52 

■ 


4.72955E-14 

2.56462E-14 

5.58290E-14 

4.69778E-14 

8.96398E-01 

8.96398E-01 

-8.96398E-01 

-8.96398E-01 

1.00094E— 00 
1.00002E- 00 

■ 


1.47465E-11 

1.47443E-11 

2.62678E-1 1 
2.C2938E-11 

2.43055E-01 

2.43055E-01 

-2.43055E-01 

-2.43055E-01 

9.989 17E-01 
1.00004E— 00 

BafPl 

1.00000E+03 

1.00000E+03 

1 .45736E-09 
1.45717E-09 

2.73856E-09 
2. 73851 E-09 

7.78425E-01 

7.78425E-01 

-7.78425E-01 

-7.78425E-01 

1.00010E-00 

1.00001E-00 

c2o7: Old Approach Timc= 6557.54. New Approach Time=5824.79 

m 

1.00000E+01 
1 .00000E+O1 

7.13873E-14 

2.69784E-14 

1 . 1 7613E-13 
4.96476E-14 

8.96398E-01 

8.96398E-01 

-8.96398E-0 1 
-8.96398E-01 

6.61955E+02 
1 .04582E+00 

■ 

1.00000E+02 

1.00000E+02 

1.4825 IE- 11 
1.47365E-11 

2.63057E-1 1 
2.62886E- 1 1 

2.43055E-01 

2.43055E-01 

-2.43055E-01 

-2.43055E-01 

7.81400E t 02 
1.08933E+00 

■ 

1.00000E+03 

1.00000E+03 

1.45761E-09 
1 .45698E-09 

2.73870E-09 

2.73832E-09 

7.78425E-01 

7.78425E-01 

-7.78425E-01 

-7.78425E-01 

6.05388E+02 

1.03607E+00 

[ c2o8: Old Approach Timo= 7424.47 , New Approach Time=6443.41 

200 

200 

1.00000E-01 
1.00000E— 01 

2.19713E-13 

2.47580E-14 

3.11883E-13 

4.59314E-14 

8.96398E-01 

8.96398E-01 

-8.96398E-01 

-8.96398E-01 

8.09897E+09 

1.23692E-05 




2.62826E-1 1 
2.63057E-11 

2.43055E-01 

2.43055E-01 

-2.43055E-01 

-2.43055E-01 

8.25499E- 09 
1.10213E— 05 


1.00000E+03 

1.00000E+03 

1.45885E-09 
1.4571 IE-09 

2.73874E-09 

2.73847E-09 

7.78425E-01 

7.78425E-01 

-7.78425E-01 

-7.78425E-01 

7.70854E— 09 
1.38970E-05 

c2o9: Old Approach Time=101 12.20 , New Approach Tiine=l 1168.74 

m 

1.00000E+01 

1.00000E+01 

8.13960E-13 

2.58682E-14 

1.04409E-12 

4.62761E-14 

8.96398E-01 

8.96398E-01 

-8.96398E-01 

-8.96398E-01 

5.13309E-17 
2.81314E- 11 


1.00000E+02 

1.00000E+02 

1.52230E-1 1 
1. 47471 E- 11 

2.65629E-1 1 
2.63039E-11 

2.43055E-01 

2.43055E-01 

-2.43055E-01 

-2.43055E-01 

5.0594 IE- 17 
4.05821E - 1 1 

20000 

20000 

1 .00000E+03 
1.00000E+03 

1. 46161 E-09 
1.45712E-09 

2.73741E-09 

2.73848E-09 

7.78425E-01 

7.78425E-01 

-7.78425E-01 

-7.78425E-01 

4.91307E- 1 7 
3.22536E- 1 1 

c2ol(): Old Approach Time— 17783.09 . New Approach Time=l 1433.52 

200 

200 

1.00000E+01 

1.00000E+01 

2.76479E-12 
2.6201 3& 14 

4.08680E-12 

4.66777E-14 

8.96398E-01 

8.96398E-01 

-8.96398E-01 
-8.96398E-0 1 

3.58499E- 25 
2.85350E* 18 

2000 

2000 

1.00000E+Q2 
1.00000E— 02 

2.05443E-1 1 
1.47421E-1 1 

2.66944E-11 

2.63037E-11 

2.43055E-01 

2.43055E-01 

-2.43055E-01 

-2.43055E-01 

2.286 14E- 25 
3.5561 IE f 18 

20000 

20000 

1.00000E-03 

1.00000E-*-03 

1.47549E-09 

1.45712E-09 

2.73275E-09 

2.73848E-09 

7.78425E-01 

7.78425E-01 

-7.78425E-01 

-7.78425E-01 

2.72428E-t 25 
3.55940E+18 

c2oll: Old Approach Time= 1803.33 . New Approach Tiino= 1869.41 

m 

1.00000E-01 

1.00000E-*-01 

9.20091E-12 

2.85327E-14 

1.09857E-11 

4.66790E-14 

8.96398E-01 

8.96398E-01 

-8.96398E-01 

-8.96398E-01 

4.76527E -■*- 32 
7.71540E 4 24 



3.32458E-11 

1.47503E-11 

4.56052E-11 
2.63034E-1 1 

2.43055E-01 

2.43055E-01 

-2.43055E-01 

-2.43055E-01 

6.81984E+32 
9.88028E * 24 


Table 2: Cost and Round-off Error Comparison of New and Old 2D Interpolation Methods With 
8 Grid Points Per Wavelength: Double Precision Computations 
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N 

T 

maxperr 

llperr 

phmax phmin 

(‘-ratio 

OLD APPROACH ABOVE. NEW APPROACH BELOW 

c2o7: Old Approach Time= 127.48 , New Approach Time=89.34 

2 

2 

1.00000E-01 

l.OOOOOE-Ol 

2.09573E-20 

2.09573E-20 

3.26076E-20 

3.26076E-20 

8.5371 IE-01 
8.5371 IE-01 

-8.5371 IE-01 
-8.5371 IE-01 

1.00000E+00 
1.00000E-I 00 

20 

20 

1.00000E-00 

1.00000E+00 

9.75947E-20 

9.75947E-20 

1.53880E-19 

1.53880E-19 

2.64616E-01 

2.64616E-01 

-2.646 16E-01 
-2.64616E-01 

1 .OOOOOE-OO 
1.00000E+00 

c2ol0: Old Approach Time=359.67. New Approach Time= 230.74 

2 

2 

1 .OOOOOE-0 1 
1.00000E-01 

1.76261E-30 

5.54668E-31 

1.51294E-30 

8.26561E-31 

8.5371 IE-01 
8.5371 IE-01 

-8.5371 IE-01 
-8.5371 IE-01 

9.99954E-01 
1.00000E- ()() 

20 

20 

1.00000E+00 

1.00000E+00 

Bii 

4.99977E-29 

8.13865E-30 

2.64616E-01 

2.64616E-01 

-2.64616E-01 

-2.64616E-01 

9.99980E-01 

1.00000E+00 

e2ol3: Old Approach Tiinc=806.21. New Approach Time= 376.07 

2 

2 

1. OOOOOE-0 1 
1. OOOOOE-0 1 

2.87010E-29 
3.38964 E-32 

4.01517E-29 

4.22501E-32 

8.5371 IE-01 
8.5371 IE-01 

-8.5371 IE-01 
-8.5371 IE-01 

1.88808E+ 16 
4.69234E+06 

20 

20 

1 .00000E+00 
1.00000E+00 

1.53766E-27 

9.39854E-32 

1.79257E-27 

1.21469E-31 

2.64616E-01 

2.64616E-01 

-2.64616E-01 

-2.64616E-01 

4.68978E— 16 
8.26702E+06 

c2ol6: Old Approach Timc=NA. New Approach Time=686.30 

2 

1 .OOOOOE-0 1 

4.3I408E-32 

3.50760E-32 

8.5371 IE-01 

-8.53711 E-01 

8. 95892E-* 28 

20 

1.00000E+00 

1.17097E-31 

1.35263E-31 

2.64616E-01 

-2.646 16E-01 

1.74818E-r29 

c2ol9: Old Approach Time— NA, New Approach Tiine=1439.74 | 

2 

1.00000E-01 

msmm 

4.17927E-32 

8.53711 E-01 

-8.53711 E-01 


20 

1.00000E+00 

1.21 719E-31 

1.51134E-31 

2.64616E-01 

-2.64616E-01 

■gfrirAMaBBill 

c2o22: Old Approach Time=NA. New Approach Time— 2195.83 

2 

20 

1.00000E-01 

1.00000E+00 

5.54668E-32 

1.21719E-31 

4.66316E-32 

1.30786E-31 

8.5371 IE-01 
2.64616E-01 

-8.5371 IE-01 
-2.6461 6E-01 

2.19468E— 76 
4.81574E-76 

c2o25: Old Approach Timc=NA. New Approach Time— 2605.29 1 

2 

20 

1 . OOOOOE-0 1 
1.00000E+00 

4.93038E-32 

1.37126E-31 

4.33575E-32 

1.55859E-31 

8.5371 IE-01 
2.64616E-01 

-8.5371 IE-01 
-2.6461 6E-01 

1.21860 r 101 

5.71654+101 

c2o28: Old Approach Time=NA, New Approach Time— 3653.88 

2 

20 

1. OOOOOE-0 1 
1.00000E+00 

5.23853E-32 

7.85779E-32 

5. 09168 E-32 
1.28541E-31 

8.5371 IE-01 
2.64616E-01 

-8.5371 IE-01 
-2.64616E-01 

2.81365 r 126 
1 .62764—127 

c2o3 

Old Approach Time=NA, New Approach Timc= 4826.52 

2 

20 

1.00000E-01 

1.00000E+00 

6.77927E-32 

7.43198E-27 

5.62612E-32 

4.69403E-27 

8.5371 IE-01 
2.64616E-01 

-8.53711 E-01 
-2.64616E-01 

3.49609+152 

3.60786+168 

c2o34 

: Old Approach Time— NA, New Approach Time— 42884.60 | 

2 

10 

12 

1. OOOOOE-0 1 
5. OOOOOE-0 1 
6. OOOOOE-0 1 

6.77927E-32 

1.25506E-16 

Infinity 

6.47594E-32 

5.60704E-17 

NaN 

8.53711 E-01 
6.01971E-01 
1.08423E-01 

-8.5371 IE-01 
-6.01971E-01 
-Infinity 

1.62494- 1 79 
1.19710 + 216 
NaN 


Table 3: Cost and Round-off Error Comparison of New and Old 2D Interpolat ion Methods With 
8 Grid Points Per Wavelength: Quadruple Precision Computations 
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7 Conclusions and Future Research 


It is now possible to develop 20(V /i order or higher wave propagation algorithms whose practical 
utility is limited only by the precision of the computer. This is made possible by a change in the 
procedure* used to implement the spatial interpolation step of the MESA schemes. \\ it h this new 
approach, the MESA schemes are far simpler to program, more efficient . and incur less round-off 
error. Automation is no longer necessary for algorithm development, but is still necessary for 
testing the MESA schemes and for applying multidimensional wall boundary conditions [4]. 

This new implementation of the MESA schemes makes it possible to adapt the method to 
the mesh instead of the common approach of adapting the mesh to t he method which adversely 
effects the CFL constraints and complicates grid generation. By adapting the method using 
estimates of the gradients from the divided difference tableau, an efficient procedure for resolving 
many different wavelength scales may be possible. In particular, problems which contain a wide 
range of wavelength scales, such as fan and jet noise generation, may benefit from adaptive 
algorithms such as these. 

In short, this new approach is an improvement to the implementation of the MESA schemes 
in every regard for methods of 15^ order or higher, and has many advantages for the 1 lower order 
methods as well. Two-point Hennitian stencils have many desirable properties [4], [7] such as: 

• They simplify boundary treatments; 

• They are more computationally efficient; 

• And. they provide higher resolution. 

In addition, they may now enable* adaptive algorithm implementations and be implemented in 
arbitrarily high accuracy using the key result from this paper, equation U). .However, a limit 
to the effective accuracy is encountered with this implementation due to the round-off errors 
incurred from subtractive cancellation as discussed in the appendix A. Some different approaches 
for reducing round-off error include: 

• Reformulate the time advance step of the MESA schemes to eliminate the factorial term 
correction discussed in section 6 which multiplies the round-off errors. 

• Modify the divided differences to central differences to avoid the divisions by grid point 
distance, h. in each column of the tableau which also magnifies round-off error [11] 

• Reformulate the MESA scheme's to propagate complex variables in time and calculate 
derivative's without subtractive cancellation [15] 

• Use divided difference pipelining [1]. 

• Use the Stirling or Bessel form of the* central differences to minimize the coefficients in 
Newton's interpolatorv equation 4. 

• Predict roundoff error using the relationship between columns of the tableau (sum of 
column j = difference of first and last terms in column j-1) [11], and t hen set to zero all 
higher order divided differences in the tableau. 

The last item limits the accuracy of the MESA scheme, but prevents contamination from 
round-off error. This also improves efficiency since it prevents unnecessary calculations which 
will only introduce error and computational overhead. The ability to detect round-off error is 
a unique advantage of this new approach that may make using very high order methods on 
available technology practical. 

It is necessary to use Hennitian stencils for very high order algorithms since traditional 
Lagrangian stencils become excessively large in the limit . This paper has provided one approach 
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for this. In addition, it has provided some of the tools necessary for exploring the capabilities 
of adaptive algorithms for adaptive 1 resolution of steep gradients (shocks). 

In the future, as very high precision, large-scale parallel systems become common. Hermitian 
divided difference implementations of the MESA schemes will be very useful since' they offer 
minimal intcrprocessor communications [4] and accuracy limited only by machine precision. For 
practical applications, this places a greater need on developing very high order wall boundary 
and radiation conditions. 
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A Sources of Round-off Error 


The benefits of very high order accuracy begin to diminish as round-off error grows. The 
predominant source of this round-off error is subtractive cancellation error which occurs when 
two numbers close in value art* subtracted from one another. This results in a loss of most of 
the significant, digits in the mantissa. 

The subtractive cancellation errors will occur in the construction of the divided difference 
tableau as the standard deviation of the values in each column decreases. Each column, num- 
bered j, in figure 3 is constructed for all i using: 

Qij = -Qi-i± ± (34) 


where for a two point Hermitian spatial interpolation, this relation reduces to dividing by x, — 
Xj-j = h. which is the distance between both grid points in figure 1. 

Each divided difference term can be proved using the generalized Hollo theorem to be equiv- 
alent to [2]: 


Q 


i'j 


j! 


-h 


< 6 , 


h 

< - 
- 2 


(35) 


That is. each divided difference j-eolumn contains data representing the j th derivative at 
some point on the interval between both points. If we knew it applied at the center (x=0) of 
the stencil, we could use it directly in the MESA scheme 1 . However, the exact location of the 
function is not known. 

From equations 34 and 35. the following is actually being calculated in each column of the 
tableau : 


f J (Si.j) = 


h 


(36) 


where £ij - 1 and are contained within the interval between both grid points, ie. xo < 

£i.j— 1 > 6 — 1 J— 1 ^ ^T ♦ 

Now let j=l and apply the mean value theorem of calculus to this equation, then it must be 
the case that : 


h 

2 


6-i.o < 6.1 < 6.o — 7^ 


(37) 


Next, let j=2. then "loosely" applying the mean value theorem again to equation 36 results in 
the relation: 

— 2 = 6-i.o < 6-1.1 < 6.2 < 6.i < 6.o = 2 (38) 

This term "loosely" is used since the denominator in equation 36 docs not change as addi- 
tional columns arc 1 processed and the mean value theorem could only be directly applied if its 
size was 6.i “ 6-i.i* However, on average, the interval will decrease for most functions since 
there will be one or more critical points (6.j) between 6-i-i and 6-i.j-J when simulating the 
oscillating functions that occur in acoustical applications. 

Therefore as j is increased, the interval in which 6 j can occur decreases. This has the effect 
of compressing all the data in column j closer together. In fact, by definition, the last column 
of the 1 tableau contains a single element representing the center of the stencil with an interval 
length of zero for the column. As the points, 6-i.j and 6.j become 1 closer in equation 36, the 
subtractive cancellation errors increase. 

For example, suppose there arc* 10 grid points per unit interval, then the first column has an 
interval length of h = 1/10. In addition, suppose the interval length is halved in each column, 
then column j will have length 6 With a 2 1 order MESA method (c2ol0) there 1 will be 
21 columns with the intervals decreasing by six orders of magnitude in the last three columns. 
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In fact, by the o th column there is very little separation between both data points (.00625 units) 
resulting in subtracting two numbers that arc* very close in value. 

This subtractive cancellation error is then magnified by the division of the small number, 
h. in equation 36. It is a simple matter to eliminate this division by using the forward, back- 
ward, or ccntered-diffcrence forms [11] of the Newton intcrpolant. However, even without this 
improvement, the results of this paper show using Hermit ian divided differences with two-point 
Hermit ian MESA schemes results in less subtractive cancellation than previous methods. 
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B Fortran90 Implementation 


This routine will calculate the spatial derivatives for the MESA scheme 
Using the data fdata(dx,dy ,x,y) = D [f [x ,y] , {x ,dx> ,{y ,dy>] 

It applies only to a c2os scheme, where the s=0,l,2, ... 

0 | 0 y=yl 


+ - y=0 


0 I o y=y0 

x=xO x=0 x-xl 

INPUT: fdata(dx ,dy ,x ,y) , x=xO,xl, y=yO,yl , dx=0,s , dy=0,s 
OUTPUT: cf(dx,dy,0) , dx=0,l, . . . , 2*s + 1 , dy=0,l, ..., 2*s+l 


program multidimensional 
implicit none 


Define variables to be used 


integer :: digits, e,i,j,k,s,n, dx, dy,r, yp, AllocateStatus 

real (kind=selected_real_kind(30) ) :: z, deltah, prodl , prod2, sum, innersum 

real (kind=selected_real_kind (30) ) :: h, x0,xl 

real (kind=selected_real_kind(30) ) , DIMENSI0N( ALLOCATABLE: : fdata 
real (kind=selected_real_kind(30) ) , DIMENSION (: , , ALLOCATABLE: : coef , Q 
real (kind=selected_real_kind(30) ) , DIMENSION , ALLOCATABLE:: cf.cfdata 
integer, parameter : : Prec30 = Selected_Real_Kind(30) 


Select the c2os MESA scheme desired. Order of scheme will be 2 * s + 1 


print *, "Enter the degree, c2os ,f 
read * , s 


Dynamically allocate memory space for array 


ALLOCATE (coef (0 : 2* (s+1) ,0 : 2*(s+l) ) > STAT = AllocateStatus) 

If (AllocateStatus /- 0) STOP ’*** Not Enough Memory ***’ 

ALLOCATE (cf (0 : 2* (s+1) , 0 : 2* (s+1) , 0 : 0 , 0 : 0) , STAT = AllocateStatus) 

If (AllocateStatus /= 0) STOP ’*** Not Enough Memory ***’ 

ALLOCATE (cfdata(0:2*(s+l) ,0:2*(s+l) ,0: 1,0:1) , STAT = AllocateStatus) 
If (AllocateStatus /= 0) STOP '*** Not Enough Memory ***> 

ALLOCATE (Q(0:2*(s+1) ,0:2*(s+l)) , STAT = AllocateStatus) 
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If (AllocateStatus /= 0) STOP '*** Not Enough Memory ***’ 

ALLOCATE (f data(0 : 2* (s+1) ,0:2* (s+1) ,0: 1 ,0 : 1) , STAT = AllocateStatus) 
If (AllocateStatus /= 0) STOP ' *** Not Enough Memory ***> 


Spacing between adjacent grid points, assume same in x and y direction 


h=dble ( 1 . _Prec30/4 . _Prec30) 


Location of right grid point as shown at top 


xl=h/2 . _Prec30 


Location of left grid point as shown at top 


x0=-h/2 . _Prec30 


Spacing between grid points xO, xl , usually same as h 


deltah = xl-xO 


! THIS SECTION NEEDS COMPUTED ONLY ONCE AND ITS RESULTS CAN BE REUSED AT EACH 
! STENCIL — SIGNIFICANT COMPUTATIONAL SAVINGS 

! The variable array, coef( , ) is assigned next and does not depend on the 
! stencil data, only upon the MESA scheme employed 

outerloopl: do dx=0, 2 * s + 1 


Compute the first set of coefficients, coef(i,dx) 


do i=dx, s 

coef(i,dx) = ( Fac (i) /Fac (i-dx) )* (-x0) ** (i-dx) 
end do 


! Compute the other set of coefficients, coef(i,dx) 

j 

csetloop: do i=s+l, 2 * s + 1 

sum = 0.0 

sumloop: do r=0, dx 

prodl=l . 0 
do e=0, dx-l-r 
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prodl = prodl * (i - s - 1 - e) 

end do 

prod2=l . 0 
do k=0 , r-1 

prod2 = prod2 * (s + 1 - k) 
end do 

sum = sum + ( Fac (dx) / (Fac (dx-r) *Fac (r) ) ) * (-xO) ** (s+l-r) * & 
( -xl) ** (i-s-l-dx+r) * prodl * prod2 

end do sumloop 
coef (i , dx) =sum 
end do csetloop 
end do outerloopl 


i 


This section is repeated for each stencil in the domain ! 


Assign known data and its derivatives from the MESA c2os 2 by 2 stencil 


do dx=0 , s 
do dy=0 , s 

f data (dx,dy, 0,0)= STENCIL DATA AT LOWER LEFT 
f data(dx , dy , 0 , 1)= STENCIL DATA AT UPPER LEFT 
f data (dx,dy, 1,0)= STENCIL DATA AT LOWER RIGHT 
f data (dx,dy, 1,1)= STENCIL DATA AT UPPER RIGHT 
end do 
end do 


Calculate the tableau by first inserting the known data into it 


lastloop: do yp=0, 1 
do dy=0,s 
do i=0 , s 


Load table in for x-direction 


do dx=0, i 

Q(i,dx)= fdata(dx,dy ,0,yp) / Fac(dx) 
Q(i+s+l,dx) = fdata(dx ,dy , 1 ,yp) / Fac(dx) 
end do 
end do 


Next perform algorithm 3.2 in Burden to construct Divided Difference Tableau 


n=(2 * (s+1)) 
do i=s+l, n-1 
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do j=i-s, i 

Q(ipj) = (Q(i,j-1) - Q(i-1, j-1))/ (deltah) 
end do 
end do 


! Evaluate the spatial derivatives at center using short form 

t 

do dx=0, 2 * s + 1 
sum=0.0 

do i=dx, 2 * s + 1 

sum = sum + Q(i,i) * coef(i,dx) 

end do 

cf data(dx ,dy ,0 ,yp)=sum 
end do 
end do 

end do lastloop 


Repeat this process for the y-direction, only using the cfdata( , , ,) 
as developed by Goodrich 


Compute Y-Direction interpolation 


Load table for Y interpolation 


dxloop: do dx=0, (2*(s+l)-l) 
do i=0, s 
do dy=0 , i 

Q(i,dy)= cf data(dx,dy ,0,0) / Fac(dy) 
Q(i+s+l,dy) = cfdata(dx,dy ,0, 1) / Fac(dy) 
end do 
end do 


! Next perform algorithm 3.2 in Burden to construct Divided Difference Tableau 

f 

n=(2 * (s+1)) 
do i=s+l, n-1 
do j=i-s, i 

Q(i,j) = (Q C i , j-1) - Q ( i-1 , j-1))/ (deltah) 
end do 
end do 


Evaluate the spatial derivatives at center using short form 
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do dy=0, 2 * s + 1 


sum=0 . 0 

do i=dy, 2 * s + 1 

sum = sum + Q(i,i) * coef(i,dy) 

end do 

cf (dx , dy , 0 , 0) =sum/ (f ac (dx) *f ac (dy) ) 
end do 

end do dxloop 


The spatial derivatives for the c2os MESA scheme at the center of the 
2 by 2 stencil are stored in the cf (dx,dy ,0,0) variables 
May need to use cf (dx,dy ,0 ,0) /(dx! dy!) depending upon how the exact 
local propagator is evaluated 


do dx=0,2*s+l 
do dy=0,2*s+l 

print *, "cf ( M ,dx, M , ,, ,dy, M ,0,0) = ,, ,cf (dx,dy,0 } 0) 
end do 
end do 


This process is repeated for the p, u, and v variables 


CONTAINS 


! Factorial Function 

! In more efficient implementations, can precalculate all the factorial 
! results and store as an array 

j 

FUNCTION Fac (N2) 

REAL (kind=selected_real_kind(30) ) :: Fac 
INTEGER, INTENT(IN) :: N2 
INTEGER : : 12 

Fac = 1.0 
DO 12 = 2, N2 

Fac = Fac * 12 
END DO 

END FUNCTION Fac 
end 
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